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The prediction of the equation of state and the phase behavior of simple fluids 
(noble gases, carbon dioxide, benzene, methane, short alkane chains) and their mix- 
tures by Monte Carlo computer simulation and analytic approximations based on 
thermodynamic perturbation theory is discussed. Molecules are described by coarse 
grained (CG) models, where either the whole molecule (carbon dioxide, benzene, 
methane) or a group of a few successive CH2 groups (in the case of alkanes) are 
lumped into an effective point particle. Interactions among these point particles 
are fitted by Lennard-Jones (LJ) potentials such that the vapor-liquid critical point 
of the fluid is reproduced in agreement with experiment; in the case of quadrupo- 
lar molecules a quadrupole-quadrupole interaction is included. These models are 
shown to provide a satisfactory description of the liquid-vapour phase diagram of 
these pure fluids. Investigations of mixtures, using the Lorentz-Berthelot (LB) com- 
bining rule, also produce satisfactory results if compared with experiment, while 
in some previous attempts (in which polar solvents were modelled without explic- 
itly taking into account quadrupolar interaction), strong violations of the LB rules 
were required. For this reason, the present investigation is a step towards predictive 
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modelling of polar mixtures at low computational cost. In many cases Monte Carlo 
simulations of such models (employing the grand-canonical ensemble together with 
reweighting techniques, successive umbrella sampling, and finite size scaling) yield 
accurate results in very good agreement with experimental data. Simulation results 
are quantitatively compared to an analytical approximation for the equation of state 
of the same model, which is computationally much more efficient, and some system- 
atic discrepancies are discussed. These very simple coarse-grained models of small 
molecules developed here should be useful e.g. for simulations of polymer solutions 
with such molecules as solvent. 

PACS numbers: 05.70.Ce, 64.70.F-, 64.75. Cd, 02.70.Tt a) Electronic mail: mognetti@uni- 
mainz.de 



I. INTRODUCTION 



It has been a longstanding challenge to predict accurately the equation of state and 
in particular the phase diagrams of fluids and fluid mixtures from atomistic models via 
computer simulation.-'^'^i^i^ Such applications have required a widespread development of 
computer simulation methodology: significant advances were possible through the inven- 
tion of Gibbs ensemble^*^"^ and configurational bias^>^»^ methodologies, grand canonical 
Monte Carlo simulations combined with histogram reweighting method3^>^»^ and finite 
size scaling^ii^'i^'i^ including field mixing,— umbrella sampling^i^ and other ex- 
panded ensemble methods.— A lot of effort has also been spent towards developing 
more and more accurate effective potentials from quantum chemistry methods (e.g. Refs. 



28 



29 



30 



31 



32 



33l ). However, for simple and industrially relevant fluids such as carbon 



dioxide^"^ it is still difficult to predict the equation of state with high accuracy, such that 
experimental data in the critical region and for temperatures ±30% around it are reproduced 
to an accuracy of a few percent.—*^ Extending such calculations to mixtures (in particu- 
lar, solutions of polymers with supercritical carbon dioxide as a solvent) is even more of a 
problem, due to the less complete knowledge of effective potentials, and due to the exten- 
sive numerical effort required. A three-dimensional parameter space involving the variables 
temperature T, pressure p and mole fraction x needs to be scanned for a binary system. 
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and the phase diagrams are typically very complicated, because vapor-liquid and fluid-fluid 
phase equilibria compete with each other.— i^i^i^ If polymers are chosen as a solute, their 
molecular weight enters as an additional fourth variable. Moreover, the coarse-grained rep- 
resentation of the solvent (e.g. carbon dioxide) and the solute have to be compatible, i.e., 
one cannot combine an atomistic description of the solvent with a much coarser represen- 
tation of a macromolecular solute. There is clearly a need to devise models that are simple 
enough to allow extensive simulation studies with an affordable effort and nevertheless accu- 
rate enough to be interesting for applications to experiment and in the context of industrial 
processing. Such validated coarse-grained models that accurately reproduce thermodynamic 
bulk properties are also a starting point for investigating the kinetics of phase separation or 
spatially inhomogeneous systems (e.g. wetting and catalysis). 

In the present work, we wish to make a step towards this goal, extending our previous 
study of a selected sample of simple pure fluids, in particular carbon dioxide^^i^ to various 
binary mixtures. We want to stress that our aim is not to reach the most accurate prediction 
of the phase diagram of a speciflc system. Indeed, motivated by the excellent results obtained 
for the pure carbon dioxide and for simple quadrupolar molecules in general,™ we want to 
investigate how this model performs for mixtures, especially solutions of various alkanes. 
In particular we will show that the new coarse grained (CG) model avoids the need for a 
big violation of the Lorentz-Berthelot (LB) combining rules (that was required in previous 
work—). This violation destroys the predictivity of the model because extensive experimental 
data for the mixture would be required to determine a parameter describing the violation 
of the LB combining rule. Due to the generality of the approach and the level of accuracy 
for the pure components,*^ the present investigation is relevant both for practical purposes 
and for a general understanding of coarse graining procedures.— i^i^ We will also present 
results of an analytical Equation of State (EOS) which (apart from some region of the phase 
diagram near critical points) is able to yield rather satisfactory predictions in agreement 
with Monte Carlo results. It is very important to note that this EOS uses the same model 
parameters as the Monte Carlo simulation. This implies that in principle we are in a position 
to attempt to predict the phase diagram of a binary mixture (which is very complex— i^i^"^) 
with comparatively small computational effort. In this view the reader should also interpret 
our choice to use LB combining rules: of course there are no reasons to believe that such 
approximations should be exact, and certainly there will be cases where more complicated 
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combining rules are preferable. However, the simple LB combining rules used here suffice 
for a wide class of systems with quite acceptable errors. 

Due to the generality of the scheme presented in this work we expect discrepancies, and 
some regions of the phase diagram might not be predicted properly. This is related to 
several limitations of the present procedure like i) the large T expansion involved in the 
building of the CG model for quadrupolar solvents^^*^ and ii) limitations related to our 
simple modeling approach like the simple potentials involved (Lennard- Jones), the neglect 
of atomistic details, and the use of the LB combining rules for which discrepancies are^ 



known to arise (see e.g. 



46l for some systems also investigated in this work). In order to 



disentangle point i) from point ii) we also present investigations of similar apolar mixtures 
for which the new CG model^*^ does not result in any improvement. The results show 
similar discrepancies f rom e xperiment as the polar phase diagrams, confirming the quality 
of the choice in Refs. 
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471 . We want to stress that in order to test the goodness of our 



CG model, the only reliable method is a Monte Carlo investigation. Indeed, without MC 
simulation it is impossible to distinguish the bias related to the approximations involved 
in the EOS from the bias involved in the CG model [point ii) above]. For instance, we 
will present results for the mixture of methane and carbon dioxide for which EOS results 
will be in better agreement with experiments than MC results: this is clearly a fortuitous 
cancellation! 

It is important to report that other interesting and significant attempts to build a sys- 
tematic description of mixture phase diagrams are present in the literature. For instance 
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49l mixtures are treated with models previously investigated in 150. For some of the 



molecules studied, these models allow for an additional parameter that can be adjusted 
and consequently a more accurate fit of experimental data is possible. On the other hand, 
there is a loss in predictivity because the full phase diagrams of the pure substances are 
required in order to determine the simulation parameters (computed in a x square fit which 
minimizes discrepancies with experiment^) plus mixture datai^"^ to determine the mixing 
parameters.— So the strategy of the present work is to deal with relatively simple mod- 
els, where (in the framework of Monte Carlo simulations) the statistical mechanics can be 
dealt with at a very good level of accuracy (e.g. long runs employing advanced Monte Carlo 
techniques are possible to minimise statistical errors and systematic errors due to finite size 
effects which are avoided by finite size scaling analysis). These models are suitable for ana- 
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lytic EOS models as well, and can serve as a starting point for the coarse grained modeling 
of polymer solutions. Of course, we do not imply that a complementary simulation strategy 
(making models as detailed as possible, to account for the packing of molecules in the liquid 
as accurately as possible, including polarizability, etc.) is not worth pursuing in its own 
right, but it is outside of the scope of the present work. 

II. COMPUTATIONAL DETAILS AND OUTLINE 

It is well establishecl^'^'^'^'^'^'^'^'^'^'^ that the most rehable approach to study the 
phase behavior of fluids is based on grand canonical Monte Carlo simulations together with 
histogram reweighting and finite size scaling techniques, especially if one wishes to include 
the critical region. In this study, we follow this approach, and amend it by successive um- 
brella sampling^ to obtain coexistence curves far from criticality. This method has the 
additional advantage that the interfacial free energy between the coexisting phases can be 
extracted as well.— "^"^"^ As we are interested in a very fast simulation code, we omit 
any potentials including effective charges, and restrict our attention to short range effective 
potentials. Three-body (nonbonded) forces are avoided as well. Electrostatic quadrupole- 
quadrupole interactions are treated as a perturbation (which is practically justifiable^), 
such that an effective angular-independent (but temperature-dependent^>^i^i^) interac- 
tion decaying proportionally to the power r~^^ of the interparticle distance r results. The 
dispersion forces are modeled by Lennard- Jones (LJ) potentials. For the sake of computa- 
tional efficiency, all potentials are cut at the distance r = = 2{2^^^)a and shifted to zero 
at Tc {(J is the range parameter of the LJ potential). When we deal with alkane chains, 
we disregard any torsional forces and bond-angle potentials and integrate a few successive 
chemical monomers into one effective monomeric unit (cf. fig. [1]). This is done in the way 
that one such unit contains three carbon-carbon bonds between successive carbon atoms, 
and we do not distinguish between interior CH2 monomers and the CH3 groups at the chain 
ends. Thus, for example, hexadecane (C16H34) is represented by a chain molecule containing 
five effective monomers (see fig. [1]).—'^ The procedure of coarsening three carbon atoms in 
a bead has been proven to be optimal in several theoretical investigations-- (see Sec. 4.3.2). 
We stress that the particular choice of coarsening three carbon units into one bead has 
nothing to do with the physical lengths of the chain (like for instance the Kuhn length), but 
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is a choice that depends more on the potentials used. Indeed as neighboring beads along 
a chain interact with a bonding potential (see Sec. IIIC for definitions) in addition to the 
Lennard Jones potential, the coarse grained model of the chain exhibits a degree of local 
stiffness, although neither bond angle or torsional potentials are included explicitely. This 
implies that the Kuhn length is longer than the diameter of our beads. 

Of course, the suitable choice of parameters is crucial for such coarse-grained models: we 
choose the strength of the quadrupole moment Q (if there is one) such that it is compatible 
with experimental data, and adjust the range a and strength e of the LJ potential such that 
the experimental critical density pc and critical temperature Tc are reproduced precisely in 
the simulation. In Sec. Ill, we will briefiy discuss the accuracy of this procedure for a variety 
of pure systems (noble gases, CO2, CH4, CeHg, short alkanes) while Sec. IV contains the 
central part of our work, in which we present a variety of results for binary mixtures. The 
additional interactions needed for the mixtures are chosen by the simple Lorentz-Berthelot 
combining rules.— 

Technical aspects of our simulations are similar to previous studies.— "^"^ Far from the 
critical point coexistence densities are computed using the successive sampling algorithm 
of Virnau and Miiller— in which high free energy barriers are overcome constraining the 
algorithm -at a certain time of the simulation- to sample configurations of a system where 
the number of particles is n or + 1. Varying n from n = to n = A^max one is able 
to reconstruct (after proper reweighting) the free energy profile F{n) at coexistence in the 
range of densities of interest. At phase coexistence, we expect a distribution F{n) with two 
peaks (corresponding to the two coexisting phases that differ in particle number) which have 
equal weight. In few very fast runs (using a small cubic box L ^ Tcm, where ctm is the 
biggest LJ length parameter of the model), invoking the equal weight rule for F{n), we are 
able to tune the chemical potential(s) to their coexistence values, with a reasonable error 
(~l-5%) which in some cases should be enough. Then, we start a second long simulation 
for a larger elongated box (to enhance the formation of the liquid-gas interface) V = 2 ■ 
with L = 9(Jm in which every window is sampled with 2-10-10'^ MC steps. Every MC 
step includes: 100 grand canonical moves in which we try to insert /delete solvent (and 
chain) particles, 1 local move in which a number of monomers equal to the total number of 
monomers are rearranged, and 10-Nchain reptation moves, where Nchain is the number of the 
chains in the box. Such a run requires on average 10 h of cpu time on 32 nodes of an IBM 
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Power4 cluster. The precision of the measured coexistence densities (for instance) is roughly 
1%. Using a spherical averaged potential allows us to speed up computations by a factor 
in comparison with the full quadrupolar model. A number of chains A^max ~ 1100 
usually allows a complete sampling of the liquid peak, while the number of solvent particles 
is typically of the same order of magnitude. We emphasize that -unlike simulations in 
the Gibbs-ensemble- in addition to the densities and compositions of the coexisting phases 
and their compressibilities, the simulation technique also provides information about the 
interface tension. At the critical point, we use the same kind of simulation described above, 
but unconstrained. (At every time the number of particles is free to fluctuate in all the 
region [0,A'max])- For more detail on the finite size analysis used we refer the reader to Sec. 
IVC. 

Even with all these simplifying approximations, establishing the phase behavior and 
thermodynamic properties of binary mixtures comprehensively still requires a lot of work 
with Monte Carlo simulations. Far away from critical points, such an effort is not needed, 
and one can try to use an analytical equation of state. We use a previously developed theory 
based on Wertheim thermodynamic perturbation theory^ (TPT). We strictly follow Ref. 



65 



66l . In particular the free energy of the system A is decomposed in a contribution due to 



a mixture of unbonded monomers (the reference system) plus a contribution due to chain 
associativity Achainr^ ^ = ^ref + ^chain- Wertheim's theory allows us to compute Achain 
perturbatively using quantities of the reference system (like pair correlation functions) and 
the known bonding potential. We use a first order perturbation theory (TPTl) which (at 
this point) reduces the problem to the computation of pair correlation functions and the free 
energy {Aj-^f) of a binary mixture of non-bonded monomers interacting with LJ potentials 
(chain-chain monomers and solvent-chain monomers) and the LJ + quadrupolar interaction 
(see Sec. IIIB) for solvent-solvent monomers. Aj-ef is computed using standard perturbation 
theory: the Ornstein-Zernike equation is solved using a Mean Spherical (MSA) closure.— 
In particular, one chooses as reference system a mixture of hard spheres with diameters 
computed using the repulsive part of the monomer-monomer potential in a Barker Henderson 
approximation,^^ while the attractive part of the potential is treated as a perturbation. A 
MSA solution is then obtained using the analytical implementation of Tang and Lu,— >^ 
in which the repulsive part of the LJ potential is fitted by a couple of Yukawa tails which 
allow to obtain an analytical result.— In our present modeling approach we need to consider 
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L J potentials plus quadrupolar interactions. This problem has been solved in |36| (see App. 
A of 36) by applying a second pair of Yukawa tails to fit the quadrupolar interaction. In 
our MSA scheme we also use a "one fluid approximation".— Our results will show how 
this simple theory is able to reproduce results in rather good agreement with MC data 
away from the critical point. On the other hand, big discrepancies occur near the critical 
points, due to the Mean Field nature of the MSA while experimental results exhibit critical 
behavior characteristic of the Ising universality class. We are aware of significant efforts to 
design proper EOS which include Ising fluctuation near the critical point.— >^ However, such 
investigations are beyond the scope of the present paper. Another popular method based 
on TPTl is known as "statistical associating fluid theory" (SAFT).— 

We want to stress that Monte Carlo simulations remain an indispensable tool in investiga- 
tions of the phase behavior of polymer solutions and mixtures. Indeed, in the present study 
the model parameters (e, a and q^) have been determined^ using the simulation critical 
points which were obtained by Monte Carlo simulation in |36|. (Any mean field approxi- 
mation has difficulties in reproducing the critical line with sufficient accuracy). Note that 
supercritical fluids are interesting and useful for practical applications, mainly due to their 
high compressibility and the concomitant large variations of density upon small changes of 
pressure, which are the origin of the breakdown of a mean field approximation like TPT. This 
means that in some very interesting regions of the phase diagram Monte Carlo simulations 
are indeed a very valuable tool. 



III. PHASE BEHAVIOR OF SELECTED PURE SYSTEMS 

When we discuss the extent to which the Lorentz-Berthelot combining rule can account for 
the phase behavior of mixtures, we need to distinguish between inaccuracies arising from an 
imperfect description of the pure components and those arising from the Lorentz-Berthelot 
rule. Therefore it is necessary to give an overview of our modeling of the pure components 
at the outset. Note that a possible additional source of errors are entropic packing effects of 
non-spherical molecules that may show up differently in a mixture of two molecules having 
different shapes rather than for a pure system, where all molecules have the same shape. 
Such effects are lost in our coarse-grained models. However, this latter criticism cannot 
be applied when we consider mixtures of noble gases, since in the framework of classical 
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statistical mechanics the description of noble gas atoms as point particles, where two such 
atoms interact with a potential depending on the absolute value of their distance only, is 
certainly appropriate. (Disregarding the case of He, quantum effects are negligible indeed at 
temperatures of interest''^). For that reason, noble gases are also included in our discussion, 
because they will bring out the possible limitations of our modeling in terms of pair-wise 
effective potentials between point-like particles most clearly. Thereafter, we shall deal with 
CO2, CeHg, CH4, and selected short alkanes. 



A. Noble Gases 



The interaction between neutral point-like particles in our work is always described by 
the Lennard- Jones (LJ) potential, 

= *^[(f )" - (f )°] ■ (1) 

Rather than working with the full LJ potential as written in Eq. ([T]), we find it compu- 
tationally more convenient and efficient to cut off this potential at r = Tc = 2''/^cr and shift 
it to zero there, such that 



Ui,{r) = Ut/{r)+AeS , f/,,(r > r,) = , (2) 

where S = 127/16384 for our choice of Tc, so that the potential is continuous everywhere. 
When we require that Eqs. ([1], [2]) yield a vapor-liquid phase diagram such that the critical 
temperature Tc coincides with the experimental critical temperature T^^^ of a particular 
system, the strength (e) of the LJ potential is fixed once T* = kBTc/e has been determined 
for the model. Likewise, requiring that the critical density pc of the model coincides with 
the experimental critical density p^^^ of that system the range (cr) of the LJ potential is 
fixed once p* = PcCr^ is known for the model. Here, T* = ksT/e and p* = pcr^ are dimen- 
sionless temperature and density, respectively. Actually, the phase diagrams of both the 
full (untruncated) L J potential and of its truncated version Eqs. [2]) have been estimated 
with high precision. '^^'-^ Fig. 11 of Ref. l36| compares these phase diagrams with each other 
and with experimental data for the noble gases Ne, Ar, Kr and Xe.''^ One can see that 
in this scaled representation the differences between the phase diagrams based on full and 
truncated LJ models are quite minor. Although noble gases are thought to be the best 
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possible experimental realization of LJ fluids, the agreement is not perfect either: while Ne 
and Ar are very close to the LJ prediction, the data for the fluid branch of Kr and Xe are 
somewhat off. This implies that even noble gases do not strictly satisfy the "law of cor- 
responding states", and hence a description in terms of classical point particles interacting 
with purely pairwise potentials of the same functional form, Uij{r) = e f{r/a), with one 
parameter for the strength (e) and another for the range (a) of the potential cannot be 
strictly true, irrespective of the form of the function f{r/a): either a more complicated form 
of the pairwise interaction, involving a third system-speciflc parameter is needed, or (what 
is usually assumed) some effect of three-body interactions^>^'^ are present. 

An even more pronounced deviation from the simple LJ model shows up, however, when 
additional quantities are analyzed, such as the vapor pressure PcoexiT) at liquid- vapor coex- 
istence and the interfacial tension 7(T) between the coexisting vapor and liquid phases of 
the ffuid (see flgs. [21 [3]). It is clear that adjusting a from p™P implies that the whole curve 
for the coexistence pressure PcoexiT) in the {p,T) plane is underestimated for both Kr and 
Xe. This is a serious drawback for the description of binary mixtures, of course, where one 
wishes to work in the {T,p, x) ensemble, x being the molar fraction of the solute. Therefore, 
we have tried an alternative, namely adjusting a such that the experimental critical pressure 
pBxp _ p^^^^{T^ is correctly reproduced. For Kr the critical temperature T^^p = 209.46 K— 
implies e = 2.8971 ■ 10"^^ J. If one uses pc = 11.0 mol/f^ to fit a one obtains a = 3.6524 A, 
while using pf^^ = 55.20 bar— instead would yield a = 3.58782 A. (For a discussion of the 
accuracy of our estimation of e and a, we refer to table I. In order to guarantee the repro- 
ducibility of our results we always present e and a with all the digits that have been used 
in our programs.) Fig. [3] shows that a somewhat better description of the vapor pressure 
Pcocx{T) is obtained over the full temperature regime from lAOK < T < T^^^. The deviation 
from the data for the surface tension a has also become smaller (fig. [2)d), but now there 
is a strong deviation between the data for the liquid branch of the coexistence curve and 
the model (fig. [2^). Similar problems are observed for Xe, where T^^p = 289. 747^^ yields 
e = 4.00747 • lO^^V, while the use of pf^ = 8.371 mo\/i yields a = 4.00053A and use of 
p^^P = 58.41 bar yields a = 3.92326A. Figs. [21 [3] show that for these noble gases the descrip- 



tion of the coexistence curve, vapor pressure at coexistence anc 
not as good as for the model of CO2 and CeHg proposed in Ref. 



surface tension is clearly 



36l . These problems carry 



over to our modelling of binary rare gas mixtures (see Sec. Ill A), as the comparison with 
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experimental data shows.— At this point we recall that as outlined in the introduction, the 
investigation of such a system has been undertaken in order to get an order of magnitude 
estimate of the errors inherent to our very simple models: the goal is not the derivation of 
a very elaborate description of noble gas mixtures. Figure [3] clearly shows that this model 
allows for a fairly good description of the mixture phase diagram (if compared to other 
mixtures presented in this work). In the present context it was not necessary to include 
more complex potentials available in the literature since long times (e.g. l8dJ8llJ82l ). 



B. Small Molecules: Methane, Carbon Dioxide, Benzene 

Methane (CH4) is also described as a point particle, and again we take Eqs. ([I]), ([2]) as 
a coarse-grained description of the interaction between methane molecules. Using T^^^ = 
130.6 K— and p^^^ = 10.1 mol/£^^ as experimental input to determine e and a, we obtain 
e = 2.63624 ■ 10^^^ J and a = 3.75792 A. Fig. H] compares the resulting model prediction for 
the coexistence curve in the temperature-density plane, the vapor pressure at coexistence 
and the surface tensions with the corresponding experimental data.— It is remarkable that 
in this case the simple potential model {Eqs. ([I]), ([2])} works better than in the case of the 
noble gas. 

For molecules such as carbon dioxide (CO2) and benzene (CeHg) the situation is more 
complicated: while CH4 is a molecule of approximately spherical shape and does not have 
a quadrupole moment, both CO2 and CgHg have quadrupole moments. Note that (at least 



36 



to a very good approximatiort^^i^) CO2 is a linear molecule while CgHg is disk-like. In 
we have shown that a very good description for both molecules is obtained when Eqs. ([1]) 
are augmented, ([2]) by a quadrupole-quadrupole interaction term. As the latter is only a 
relatively small perturbation of the Lennard- Jones-type interaction, it suffices to treat the 
(angular-dependent) quadrupolar interactions via thermodynamic perturbation theory. To 
leading order this yields the following effective potential^'^'^ 

Here, Q is the strength of the quadrupole moment of the considered molecule. Note that 
the interaction is isotropic and inversely proportional to temperature. We also cut off this 
part of the interaction at the same radius Vc as the LJ interaction, and shift it to zero at Tc 
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as well, which yields the following total pairwise interaction for these molecules 



^^^^^^ ^ . 4e[(a/r,,)i2 - (a/r.,)^ - M^/r.,r + S],r<r^ 

, r > Tc 

where 

- .7q 
^ ~ 16384 + 5 256 

and q is the reduced quadrupolar interaction parameter, 

q = Qy[ea''kBT]=q,TjT , qc = q{T,) . (6) 

Note that Eq. is given in CGS units; in SI units, there would be an additional factor 
(4vreo)-2. 

Using Eqs. Q-Q, one can fix e and a such that critical temperature T^^^ and density 
p^^'P are reproduced. (For Q, the experimental value is taken as a first guess). As discussed 



in Ref . 



36 



this leads to a self-consistency problem, since Eq. ([6]) must hold together with 



e(ge) = knTr/T:iq.), AQc) = [^^^1^] , (7) 

Pc I^A 



where MmoI is the molar mass of the molecule and Na is Avogadro's number. This 



problem was solved in |36| by determining the functions T*{qc)/T*{0), and p*{qc)/ pciS^) by 
extensive Monte Carlo simulations for a broad range of values for q^.. It turns out that for 
CO2 the experimental value Q = 4.3 ± 0.2DA yields 

qc = 0.387, e = 3.491 x 10"^^ J, a = 3.785 A , (8) 
while for the case of benzene the value Q = 12DA would imply 



qc = 0.247, e = 6.910 x 10"^^ J, a = 5.241 A . (9) 

The corresponding results for the vapor-liquid coexistence curves in the (T, p) and {p, T) 
planes as well as the temperature dependence of the interfacial tension for both CO2 



and CgHg were already presented in |36| and shown to give a rather good agreement with 
experiments.— 
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Of course, the disregard of the angular dependence of the quadrupolar part of the in- 
teractions is a matter of concern. This point was investigated by us in 37(, where detailed 
comparisons of Monte Carlo results for the full angular-dependent quadrupole-quadrupole 
interaction and the isotropic approximation {Eqs. (l3])-(j6l)} were performed for the case of 
CO2. It was shown^ that the model with LJ + full quadrupolar interactions (which is still 
a crude coarse-grained model, in comparison with all-atom models including partial charges 
etc.) does not provide a better account of the experimental data than the spherically aver- 
aged one. 

Another point of concern is the possible sensitivity of the results of such models to 
the precise value of qc. Note that q is proportional to {Eq. ([6])}. Consequently, a 
small experimental error in Q is magnified considerably. There may also be systematic 
effects since Q is often determined in the dilute gas phase. Here, we are interested in using 
densities around the critical density, and Q could be slightly renormalized there. Packing 
effects should also be taken into account. Indeed, CO2 is not a spherical molecule, and 
at high density a local orientational order could arise. This packing could enhance some 
favorable angular correlations that give rise to a higher effective quadrupolar moment. One 
can argue that high temperature perturbative theory {see Eq. (|3])} may not be very accurate 
and higher order terms could be important: in fact our previous investigation in which a 
full (angular dependent) quadrupolar interaction was considered,— proves that this is not 
the case. In addition, one may argue that the model of Eqs. (I3])-(I6]) is an effective model, 
intended for a good representation of equation of state data, particular for vapor-liquid 
equilibria (VLE). Therefore, qc should be treated as an effective parameter which can be 
used to optimize the description of such VLE data. In this spirit, we have also tried different 
choices of qc and found that a slightly better description of CO2 is obtained 



This choice was already included in our previous work.— 1^ For benzene, a very good 
agreement with experiments can be achieved for 



Fig. [5] presents the coexistence curve of benzene in the p — T and T — p planes as well 



= 0.47 , e = 3.349 ■ lO'^^ J, a = 3.803 A. 



(10) 



= 0.38 , e = 6.472 ■ 10"^^ J, a = 5.284 A. 



(11) 
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as the interfacial tension. Results based on Eq. (fTT!) are compared with results based on 
the previous choice {Eq. Q}^^ and with experimental data.— The description of the ex- 
perimental data is clearly remarkable over a wide range of temperatures. It turns out (see 
below) that these "optimized" choices of parameters {Eqs. ( |T0|) . ( ITTi) } also yield a much 
better description when we consider mixing behavior (e.g. CeHg + CH4). 

C. Short Alkanes 

In this section we briefly discuss the extension of our methodology to systems such as 
propane (CsHg), pentane (C5H12) and hexadecane (C16II34). These short alkanes are just 
treated as test systems for our methodology and will be used in Sec. Ill as components 
in binary mixtures. Our methodology can be used, in principle, for any alkanes, provided 
information on the vapor-liquid critical point (T^^P,p^^P) is available. (Unfortunately, this 
is not the case for much longer chains). 

As it was already emphasized (fig. [1]) we do not attempt an all-atom description of alkanes. 
We also do not use an united atom model where CH2 (or CH3) groups are described as one 
spherical pseudo-atom.— 1^ Such a model requires torsional and bond angle potentials and 
is still rather demanding to simulate. As indicated in fig. [H we reduce the description to a 
coarse-grained bead-spring model, where a small number of successive CH2 or CH3 groups 
are combined into a single effective monomeric unit. For C16H34 we choose 5 effective units, 
so each unit contains about 3 C-C bonds. For pentane and hexane we choose a dimer (but 
the effective LJ parameters e and a are different, of course). Such a model is perhaps most 
questionable in the case of C3H8, which we treat as a single effective unit (i.e., such molecules 
are treated like almost spherically symmetric molecules such as methane). 

We keep the (truncated and shifted) LJ potential {Eqs. ([T]), ([2])} between all pairs of 
effective units, bonded and non-bonded ones. In addition we use the well-known FENE 
potential for the bonded ones^ 

t^FENE(r) = -33.75e ln[l - [r/l.baf] (12) 

We note that in Eq. ( fT2|) e and a are the same parameters as in the LJ potential between 
the monomers. The parameters of the FENE potential have been chosen to prevent the 
crossing of macromolecules in the course of their motion. We note that this choice does not 
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reproduce the characteristic ratio of alkanes accurately. This means that the FENE potential 
is fully constrained, and the model remains a two parameter model with parameters chosen 
to match the critical temperature and density. 

On this coarse-grained level both torsional potentials and bond-angle potentials between 
effective beads are ignored. Hence, it is worthwhile to test whether such crude models are still 
able to reproduce the phase diagram and other thermodynamic properties of the real system 
correctly. Thus, figs. [MHl show results for the phase diagrams of several members of the alkane 
series (including CsHg, C5H12 and C16H34) in the T-p plane, as well as the corresponding 
coexistence pressures and interfacial tensions between the coexisting vapor and liquid phases. 
The agreement between the model results and the corresponding experimental data^ is 
remarkable, again, although it is not as convincing as for methane (which we have included 
for comparison). In particular, for C5H12 deviations clearly occur. Table I collects the 
experimental critical temperatures, densities, and pressures,— as well as our choices for e 
and a for the materials studied, and the prediction for the critical pressure that results from 
our model. 

In all cases the critical pressure is predicted with an accuracy of a few percent, and a 
glance on Fig. [7] shows that the slope of the vapor pressure versus temperature curve is close 
to the slope derived from experiments, too. For temperatures away from the critical region 
(say, 20% below T^), deviations between experiment and the model predictions become 
visible, both in the coexistence curve, coexistence pressure, and interface tension (Fig. [HI), 
in particular for propane and pentane. Of course, the accuracy of the modeling could 
be enhanced by allowing for additional adjustable parameters like in many models in the 
literature, e.g. by introducing a bond-angle potential, or more interaction sites (see e.g. l45l). 
Then, quantities such as the acentric factor (referring to the shape of the coexistence curve 
30% below Tc*^) can presumably be fitted nicely. However, the simplicity of the coarse- 
grained model is lost. Experience with such somewhat more complicated models shows that 
these models still require correction parameters to the LB combining rules that deviate 



from unity by about 10% (see e.g. [49|). Without these additional parameters (note that 



it is not at all straightforward to find optimal values for these parameters) the gain in 
accuracy that such models yield for the description of mixtures is rather modest. Note that 
an important motivation for the present work is to develop simple models suitable for the 
simulation of polymer solutions (the case of hexadecane in CO2 being just a prototype case). 



16 



We are not focusing on pushing the accuracy of modehng of pure short alkanes to its hmit. 

IV. PHASE BEHAVIOR OF SELECTED BINARY MIXTURES 

Extending our treatment to binary systems (A,B) one wishes to describe the interactions 
between unhke particles by a potential of the same functional form as it is used for the 
interactions between particles of the same type, i.e. the Lennard- Jones potential in our case. 
The simplest choice, most often used in the literature, is the Lorentz-Berthelot combining 
rule^ 

0"AB = (o^AA + CrBB)/2, eAB = \/^AA^BB (13) 

As is well-known, there is really no convincing derivation of Eq. (fT3l) . so there is no 
reason to believe that Eq. f |T3l) is exact. At best it is a practically useful approximation. As a 
matter of fact, several alternatives to Eq. ( fT3l) have been proposed in the literature.— '^'^'^'^ 
Although it has been demonstrated that there are some cases where some of these alternative 
combining rules work better, in general none of these alternative combining rules has a really 
clear advantage.-^ Since we wish to explore a very simple and general approach, we do not 
implement any alternatives to the simple Lorentz-Berthelot rule in our paper, even when 
one has to pay the price of sacrificing a small improvement in the accuracy of our modeling. 
We also note that the Lorentz-Berthelot rule works very well for the prediction of virial 
coefficients for the mixture of Argon plus CO2, a mixture of an apolar and a quadrupolar 
fluid.^4 We want to stress that proceeding in such a way no experimental input from the 
mixture phase diagram is required for testing a full predictive model for the mixture. This 
also holds for the TPTl computations which require only e and a that can be obtained using 
Monte Carlo results of the pure component critical line.— Coexistence densities and pressure 
have been computed as in pure component systems.— On the other hand, the computation 
of the critical points is more complicated. Indeed, in a binary mixtures close to criticality 
the proper identification of the order parameter is a subtle problem.— In principle, complete 
scalingS^i^ii^ in the case of binary mixtures implies that three scaling fields occur, which are 
linear combinations of four independent intensive variables: the deviations of two chemical 
potentials, temperature, and pressure from their values at the critical point. Consequently, 
the order parameter density becomes a function of the appropriate conjugate variable, and 
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the relevant physical densities (particle number densities, entropy density) become nonlinear 
functions of the proper scaling fields.— Since this formalism is somewhat cumbersome for 
the case of compressible binary fluid mixtures, we simplify the problem by applying "field- 
mixing" -procedures analogous to the method of Wilding^^'^'^ which is rather successful 
for most one-component fluids. Details on this procedure are reported in the appendix 1X1 
presenting the analysis done for a critical point of the Krypton Xenon mixture. In order 
to estimate systematic errors of this procedure, in appendix |X] we also present results with 
a full finite size analysis with cumulants crossingi^ for a highly asymmetric mixture like 
carbon dioxide in hexadecane. 

A. Mixtures of Small Apolar Molecules 

As a first example of apolar mixture we present results for kripton plus xenon. As it has 
been discussed in Sec. II, the noble gases already exhibit rather large deviations between the 
experimental data and the model calculations based on the Lennard- Jones potential. Thus, 
it is interesting to see whether these problems get even worse when mixtures are considered. 
The resulting critical line in the {p, T) plane is shown in fig. [3] for both choices of e and cr as 
discussed in Sec. II A. If we fit pc and Tc for the pure systems, the predicted critical points 
for the mixture deviate from the experimental curve about as much as for the pure systems. 
If we adjust e, a such that Pc,Tc is reproduced, the datar^ for the two mixed systems are 
almost perfectly reproduced. The variation of the critical concentration with temperature 
is also rather well reproduced (fig. [9]) by both models where pc and or and are fitted 
to experimental values. 

As a second case we consider now methane in butane. In Sees. II B, C we showed 
that the simple LJ model gives a fairly accurate account of the equation of state of both 
CH4 and CsHg. Therefore, it is natural to consider a mixture of those two molecules as 
a next step. Of course, a comprehensive study of the phase behavior of such mixtures in 
the space of all three variables {T,p,x) is a nontrivial effort. Therefore we limit ourselves 
to consider only isothermal slices through the phase diagram, following a standard practice 
in the literature.—*^ As an example, fig. [10] shows two such slices at T = 327K (a) and 
T = 277 K (b), and compares experimental data^^ with selected Monte Carlo data and 
results from our implementation of the TPTl-MSA (which is described in appendix B of^^). 
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We emphasize that the various parameters characterizing the interactions among the various 
molecules are those obtained from Monte Carlo simulations of the pure materials (Sec. II), 
together with the Lorentz-Berthelot rule. These parameters also serve as input for TPTl- 
MSA: there are no additional parameters that enter the latter approach. Thus we present 
comparisons between experiments, simulations and theory in which no adjustable parameters 
for the mixture have been used. 

It should be noted that both chosen temperatures in fig. [10] fall below the critical tem- 
perature of CaHg but exceed the critical temperature of CH4. Therefore, the characteristic 
bubble-shaped liquid-vapor coexistence curve results, starting out at the ordinate axis at 
the vapor-liquid coexistence point of pure CsHg, but not extending to CH4 concentrations 
close to X = 1. The critical point occurs at the maximum of this closed loop. (The liquid 
phase is located on the upper part of the loop to the left of the critical point, the remaining 
part of the curve describes the vapor). For T = 327K and x < 0.35 both experiment, TPTl 
and Monte Carlo agree nicely. For larger x, however, a systematic discrepancy between 
Monte Carlo data and experiment shows up. The TPTl-MSA approximation overestimates 
the critical pressure substantially. This problem already occurs in the pure systems, as is 
well-known, and is an inevitable consequence of simple mean-field-like approximations.—'^'^ 
Fig. m shows that the critical temperature and pressure of pure CH4 are both overestimated. 
The same holds for pure CsHg, and the whole line of critical points Tc{x) that connects 
Tc(0) and Tc(l) when we would project them into the {p,T) plane as we did for the Kr-Xe 
mixture (fig. [3]). As in the latter case, the mixture of CH4 and CsHg has a simple "type 
I" phase diagram in the classification scheme of fluid binary mixtures^'^'^ (type 1^ in the 
modern classification'^^). As a consequence, we expect that TPTl-MSA predicts too large 
vapor-liquid coexistence loops in the {p, x) plane at all temperatures that are supercritical 
for CII4 but subcritical for CsHg. 

A more disturbing discrepancy seems to occur between the datai^ and the theoretical 
results at the lower temperature (T = 277 K), where at small x the vapor pressure at coexis- 
tence falls slightly but systematically below the experimental data. For molar concentrations 
well below criticality, Monte Carlo results and TPTl-MSA agree very well, and our numer- 
ical procedures are accurate for our model. Hence, assuming that the experimental data 
are accurate enough so that the discrepancy is meaningful, this result indicates that some 
limitations of our model become apparent. This is not really a surprise, of course, because 
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in the data for pure propane at this temperature discrepancies of the order of a few percent 
do occur as well (figs. [MH])- 

As a third case we now consider the mixture of CH4, and C5H12, because for pentane 
slightly larger deviations between the predicted and observed coexistence vapor pressure do 
occur over a much broader temperature range (fig. [7j). Indeed, the corresponding isothermal 
slices through the phase diagram of that mixture (fig. [TT]) . which still is a type-I phase dia- 
gram, show that slight but systematic discrepancies are now seen at the higher temperature 
as well. At the low temperature, the phase diagram can only be reproduced in a rather 
qualitative manner. Note, however, that T = 237i^ is less than 50% of the critical temper- 
ature of pentane, where the effective interactions of pentane were adjusted: of course, the 
coarse-grained modelling used in our work should not be pushed to too low temperatures. 
Keeping this limitation in mind, we conclude that a rather satisfactory description of mixing 
behavior of these systems is in fact reached by our models. Hoping for perfect agreement 
would have been premature, in view of the simplicity of our models. But the phase diagram 
predictions should allow a useful first orientation at temperatures not too far below of the 
higher critical temperature of the components in such a binary mixture. 

B. Mixtures of small molecules, one of which has a quadrupole moment 

We begin with a mixture of CH4 and CO2, because for both pure molecules a particularly 
accurate description of the equation of state was obtained (see Sec. II). Again we note 
that the CII4 + CO2 system belongs to the category of "type I" phase diagram in the 
classification scheme of Scott and van Konynenburg^i^i^ (l'^ in the modern classification:^) 
and the temperature regime of interest for our modeling is the regime in between the critical 
temperatures of the two constituents of this mixture. Note that Eq. (JT3i) only applies to the 
LJ part of the interactions of CO2, since CH4 has no quadrupole moment. 

In fig. [12] we present isothermal slices through the phase diagram in the space of variables 
{T,p,x). If one uses TPTl-MSA the model for CO2 based on Eq. (ITOl) can describe the 
mixing behavior with CH4 very accurately at molar concentrations x of CH4 and pressures 
that are not close to criticality. As emphasized above, mean-field theories such as TPTl- 
MSA are not expected to be accurate near critical points. Hence, the discrepancy that TPTl- 
MSA predicts a too large loop inside of which two-phase coexistence occurs, is inevitable and 
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expected. But for the model Eq. (flOi) the part of the loop at not too large x is significantly 
more accurate (full curves) than a simple LJ model for CO2 would be (broken curves). As 
expected, at low temperatures (such as T = 230/^) the quadrupolar model for CO2 {Eq. (IHl)} 
also starts to show slight but systematic deviations from the experiment at the vapor branch 
of the vapor-liquid coexistence curve. This is similar to our finding for the apolar mixtures 
(Sec. Ill B). 

In order to verify that the good agreement between experiment and theory for the 
quadrupolar model of CO2 in the CH4+ CO2 mixture is not just fortuitous, we show in 
fig. [13] corresponding results for the mixture of benzene (CeHg) and methane (CH4). This is 
a more stringent test, since the critical temperatures of the two constituents are rather far 
apart from each other (cf. figs. HI [5]). Nevertheless, the conclusions are the same as in the 
case of CH4+CO2: using interaction parameters that were optimized for the pure systems, 
namely those of Eq. (fTT!) in the case of CgHg, and adjusting them to Monte Carlo results as 
described in Sec. II, we can proceed to the description of the mixture datar^ and estimate 
the missing mixed interaction parameters from the Lorentz-Berthelot rule, Eq. ( fT3l) . The 
use of these interaction parameters in a simple and fast analytical theory for the EOS such 
as TPTl-MSA then provides a satisfactory description of the phase behavior of the mixture, 
apart from the vicinity of critical points (this drawback can be rectified by carrying out MC 
work for the mixture as well, of course) and for not too low temperatures. (For tempera- 
tures of the order of 50% of the critical temperature of the constituent with the higher 
Tc systematic deviations start to appear rather generally.) 

The last example of this section deals with a slightly more complicated case, namely the 
CO2 + C5II12 system (fig. [14]): while CO2 is still represented as a point particle with a 
quadrupole moment, as in the previous examples, the other partner of this mixture (C5H12) 
should not be coarse-grained into a point particle any more, but rather needs to be rep- 
resented as a dimer (i.e., a dumbbell-like effective molecule). In this case the TPTl-MSA 
theory predicts unmixing over a far too large range of molar CO2 concentrations, and the 
improvement provided by the inclusion of the quadrupolar moment at small x is only qual- 
itative, but not quantitative. On the other hand, the Monte Carlo results for this model 
are in rather good agreement with the corresponding experimental data.^''^ Since MC and 
TPTl-MSA are using precisely the same interaction parameters, we conclude that for this 
particular case TPTl-MSA is somewhat inaccurate for the vapor branch of the mixture. 
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far away from criticality. A related discrepancy was already noted for the CH4 + C5H12 
system at T = "ilSK (fig. [TTb). Perhaps this indicates that TPTl-MSA does not capture 
the statistical mechanics of flexible dimers well enough. 



C. Polymer solutions: The CO2+C16H34 system revisited 

Virnau et dX.—'^^^ already attempted to model this system, describing CO2 as a point 
particle with no quadrupole moment. They found that using the Lorentz-Berthelot rule 
{Eq. ( fT3i) } the phase diagram predicted by the model belongs to type I, while experiments 
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where 



suggeslr^i^ that this system belongs to the type III class according to 

1^ means that the critical line emanating from the pure component critical point of the 
hexadecane goes to high pressure regions without joining the solvent critical point like in 
diagrams starting with 1^) . Virnau et al.— proposed that one can improve the description 
by using an empirical factor ^ to modify Eq. ( fT3|) . assuming that sab = C,\^^aa^bb instead 
of eAB = \/^AA(-BB- In the literature, the value of ^ depends on the specific mixture and 
typically is written in the form ^ab = 1 — ^ab-, with kAB > 0. Of course, there is not really 
a theoretical justification for doing so, and ^ simply plays the role of a fitting parameter. 
By trial and error it was found that ^ = 0.886 provides a description compatible with the 
experimental data. 

In the present subsection of our paper, we show that the main source of the problems 
encountered in|42| was the neglect of the quadrupole moment. Thus, we have repeated the 
study of the CO2+C16H34 system, insisting on the Lorentz-Berthelot rule, Eq. f|T3l) . but 
using Eq. (ITO!) as an improved model for CO2, as in the previous subsection. Again, the 
Lorentz-Berthelot rule is only applied to the Lennard- Jones part of the interactions, since 
C16H34 does not have a quadrupole moment. 

Following the strategy of the previous subsection, we have computed an isothermal slice 
through the phase diagram at T = 4867^^, where data from the previous simulation^ were 
available both for ^ = 1 and for ^ = 0.886. Indeed it is found that the data of the present 
model = 1, but optimized quadrupolar interaction = 0.47 for pure CO2) are well 
compatible with the experimental data^°^ and almost fall on top of the results of the previous 
calculation with qc = and ^ = 0.886.^- 

Of course, we have already seen in the previous subsections, that often a very good 
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agreement between our description based on a very simplified model occurs at high enough 
temperatures. In order to test, to what extent this problem arises for the present system, 



we have followed the strategy of [4^ to compute the full critical line Tc{x) , pc{x) for the full 
range of molar concentrations x of CO2. Fig. [T6] shows the resulting projection into the p*, T* 
plane. (Here p, T are given in LJ units, with the L J parameters of the effective monomers 
used to rescale the variables). One sees that the simulations with nonzero quadrupole 
moment included in the figure are close to those for ^ = 0.9, Qc = 0, for T* < 1.3. As 
a consequence, the model that we have developed for CO2, Eq. (fTOl) . is still not able to 
yield the correct phase diagram topology. (For ^ = 0.9, transition type IV was observed 



in Ref. 
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as opposed to type III which was observed experimentally.) For T* < 0.8 
the model does not yet describe the properties of hexadecane + carbon dioxide mixtures 
accurately, although for T* > 1.3 (T > 545K) the properties of the system are predicted 
rather satisfactorily. Of course, this result is not unexpected. For T < O.bT^"^^ ^ 360K the 
model based on fitting the critical parameters of hexadecane to fix its interaction parameters 
starts to become inaccurate. On the other hand, the proper prediction of the phase diagram 
type is a very stringent test. Indeed, variation of the interaction parameters by a few percent 



could drastically change the type of the phase diagram 
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V. CONCLUSIONS 



In this paper, we have studied the phase diagrams of a variety of fluid binary mixtures, 
with particular emphasis on mixtures of alkanes in supercritical carbon dioxide and benzene. 
In order to better understand the performance of our modeling for these systems, we have 
also investigated mixtures with apolar solvents including noble gases and methane. We have 
investigated the accuracy of the use of the Lorentz-Berthelot rules for describing the mixing 
behavior, based on interaction parameters for the pure systems that are tuned such that 
the critical point (critical temperature, critical density or pressure) of the pure systems are 
well reproduced. Using a simple Lennard- Jones model for interaction parameters of pure 
apolar fluids, Monte Carlo calculations in the grand-canonical ensemble, analyzed by appro- 
priate finite size scaling methods, readily yield the desired accuracy for this procedure. For 
the polar molecules we use a spherically averaged point-like quadrupolar interaction,—'^'^ 
which was shown to produce very good phase diagrams, also if compared to more realistic 
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atomistic models. Our model takes as experimental input the critical temperatures and 
densities of the pure components (like in previous coarse grained schemes^) plus the exper- 
imental quadrupolar moments. For pure CO2 and CeHg this choice leads to a significant 
improvement in comparison with a simple LJ model without explicitly accounting for the 
polar interactions. In Ref. l36| and Fig. [5], as a second option, we have treated the quadrupole 
moment as an effective parameter in an attempt to optimize agreement with experiments. 
We tune this parameter such that the liquid branch of the vapor-liquid coexistence curve 
of pure CO2 or pure CgHg is optimally represented. In the case of benzene, for which the 
optimization procedure seems to work very well, the agreement with the coexistence pres- 
sure is also improved. (This is not the case of CO2 which is however better described than 
benezene if the experimental values for the quadrupole moments are used.) The physical 
reason for this requirement to work with an effective quadrupole moment is presumably that 
actual molecules are not point-like particles, of course: CO2 is a rather elongated molecule, 
while CeHg is disk-like. So packing effects should occur, i.e. local orientational correlations, 
which are underestimated by the quadrupolar interaction. Thus, it is gratifying to note that 
a remarkable improvement of accuracy in the prediction of the phase behavior of mixtures 
is achieved if this effective quadrupole moment is used. 

These energy parameters, which we fixed from the description of the pure systems, to- 
gether with the Lorentz-Berthelot rules, allow us to predict phase diagrams of mixtures, 
with no ambiguity whatsoever, since no further adjustable parameters occur. Two methods 
of prediction are used: (i) Monte Carlo simulations (ii) TPTl-MSA calculations. The Monte 
Carlo approach has the substantial advantage that it is also accurate near critical points of 
the mixture. In principle, we obtain the exact statistical mechanics of the model system. 
Any discrepancy between experiment and prediction is entirely due to a shortcoming of the 
(simplified) model. The TPTl-MSA approach has the merit that relatively little computa- 
tional effort is necessary to implement it. However, it clearly involves various approximations 
and hence the interpretation of discrepancies between TPTl-MSA and experiment is not so 
clear - part of them being due to inadequacies of the model, part of them stem from inaccu- 
rate approximations. For instance, TPTl-MSA, like all mean-field theories, overestimates 
the critical temperature and pressure, so the isothermal slices through the phase diagram of 
the mixture always involve two-phase regions which are too large. 

We note that fiuids like CO2 have an important application as supercritical solvents. If 
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one aims at describing the behavior in the critical region of the pure solvent and of the 
mixtures correctly, Monte Carlo methods have a clear advantage. Now, one could try to 
readjust parameters in the TPTl-MSA approach to improve agreement with experiment (like 



done for instance in 
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where e and a for the EOS have been rescaled in comparison to the 
model used in MC simulations in order to properly reproduce the critical points of the pure 
compounds), but this would be just an attempt to provide a partial cancellation of errors, 
and in other regions of the phase diagram the description would necessarily get worse.— 
Since we feel that relatively little physical insight is gained by such fitting procedures, they 
have not been implemented in our paper. Our overall conclusion is that in the framework of 
the modeling as defined above the Lorentz-Berthelot rules work very well, in the sense that 
an ad-hoc change of mixed binary interactions by at most a few percent (typically one or two 
percent) would lead to almost perfect agreement with experiment. As a piece of evidence for 
this claim, we note that in the study of the CO2+C16H34 system by Virnau et al.,— where 
the CO2 molecule was modeled as a point particle with LJ interactions with no account of 
the quadrupole moment, a correction factor ^ = 0.886 to the Lorentz-Berthelot rule was 
required to produce good agreement with experiment. However, the present model (with 
a quadrupolar interaction and no correction factor) yields results that are almost identical 
to those of Virnau et al.-^ when ^ = 0.900 is chosen. As a consequence, we conclude that 
in the present model a correction factor ^ ^ 0.985 would suffice to reproduce the results. 
Noting that the Lorentz-Berthelot rule assumes that the mixed correlation functions in the 
fluid described atomistically behave in the same way as in the coarse-grained descriptions, 
deviations of such a fitting parameter ^ from unity in the range from 1 to 2% are no surprise 
at all. 

Thus we feel that the present level of accuracy cannot easily be improved in the frame- 
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complicated models for fiuids are discussed in the literature (see 
timizing parameters in those models such that the critical properties of the pure systems 
and their vapor-liquid coexistence curves are very well reproduced, might be an alternative 
starting point to test the validity of combining rules.— 1^1^ However, even for our very sim- 
ple model the Monte Carlo runs require substantial computer resources, and hence we have 
not attempted to generalize our approach to other models. The strength of our method, if 
compared to more detailed models with a lot of parameters the optimization of which would 
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require massive computation, is its generality. It is possible to have the potentials for a given 
mixture, without any extra computational efforts from the results of the pure components.— 
It is also important to mention that the present work validates the use of spherical averaged 
quadrupolar potential.— i^"^ 

Finally, it is important to observe that our way to coarse grain solvent molecules into 
single beads has the advantage, with respect to atomistic models like multi center Lennard 
Jones, to be accessible to advanced equation of state machineries. In this paper we have 
shown how, with rather small efforts, significant results can be obtained. We are also aware 
of the fact that several improvements could be done (e.g. using some integral equation 
scheme which should improve the MSA solution near the critical point). However, most of 
the advanced methods in equation of state modeling, apply only for reference systems that 
are mixture of monomers (i.e. beads with point-like interactions), the associating part being 
taken into account by TPTl. On the other hand TPTl gives a reasonable description only 
if in the "associated molecule" diameters of the beads do not overlap. This is of course the 
case in our CG model for alkanes in which the experimental distance between three carbon 
units (d=4.59 A) is bigger than the typical a used (cr ~4 A). (This is another reason the 
FENE potential uses the same simulation parameters of the LJ interaction.) The condition 
d> cr guarantees that the reference system (a mixture of monomers) is a good starting 
point for a perturbation theory. On the other hand if d< a association is too strong and 
cannot be properly taken into account by TPTl (i.e. the monomer reference system is not 
the adequate starting point for a perturbation expansion). In models which describe simple 
solvent molecules with several interacting points (see e.g. Tab. 1 of |50| for typical parameters 
of two center LJ) we have d< a: this implies that these models cannot be investigated with 
associating theories. In conclusion, our modeling approach might enable the application of 
modern EOS which is another important motivation of our work. 

We do feel that the approach based on the present models is able to make nontrivial 
and practically useful predictions for a large class of systems. Hopefully more experimental 
data on mixtures will become available to allow for more stringent tests. 
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APPENDIX A: DETERMINATION OF THE CRITICAL POINT 

In this appendix we describe how the critical points for the mixtures (e.g. figs. [3]fT6|) have 
been obtained. As mentioned in the sec. [IV]in asymmetric binary mixtures the determination 
of the order parameter is difficult^ and a complete scaling^^i^ZiS^ requires a lot of work which 
has to be repeated for every critical point. For this reason in the present work the method 
of Wildingi^i^Siiii has been used. Taking as an example the kripton xenon mixture (but the 
same procedure has been done for all the other mixtures), we take the order parameter M 
as a linear combination of the particle numbers of Kr atoms {Nk) and of Xe atoms {N^) and 
of the total potential energy Etot 

M = N^ + xiNk + X2Etot , (Al) 

where the parameters Xi,X2 are determined by the following iterative procedure. First, the 
chemical potential /x* (in LJ units) is tuned to get vapor-liquid phase coexistence (i.e., the 
distribution P{M) satisfies the equal area rule). Normally, due to the lacking of the particle- 
hole symmetry, the two peaks will not be symmetric such as the two peaks of the universal 
Ising model distribution at criticality. Thus, as a second step the other chemical potential 
is also tuned (and the step 1 is repeated), so that one gets somewhat closer to the critical 
point of the system. Still, the universal shape of the distribution is not yet obtained. The 
third step consists in a variation of Xi (and repeating steps 1 and 2) such that the two peaks 
of the distribution become as similar to each other as possible. The fourth step amounts to 
a variation of X2 (again repeating steps 1 and 2). In this way (for the investigated cases) it is 
possible to obtain a final histogram of the order parameter M that reproduces the universal 
Ising shape at criticality almost exactly (fig. fTTl) . 

On the other hand the previous approach is not totally correct because it neglects finite 
size corrections for the critical parameters but simply "supposes" that the simulation box is 
large enough. In order to elucidate this point in this appendix we also report the full finite 
size analysis^^ with crossing cumulants for several simulation boxes. We do this investigation 
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for the polymer solution studied in this paper (CO2+C16H34) which should be more sensible 
to mixing effects being a highly asymmetric mixture in which the coupling between two 
order parameters (total density and relative concentration, respectively) may be more of a 
problem rather than for some noble gas or small molecule mixtures. In fig. [18] we report our 
finite size analysis for the critical point at T=486 K. From extensive figfJ-pVT simulations (/i^ 
being the chemical potential of the solvent, CO2 in this case, and /ip the chemical potential 
of the polymer C16H34) histograms have been generated for the total energy Etot, the number 
of solvent particles A^^ and the number of polymers, Np. The probability distribution for the 
order parameter M = Np + xiNs + X2Etot is computed, using xi = 0.08 as an initial guess (it 
turns out that the final results depend on the parameter X2 so weakly, that one may choose 
X2 = here; on the other hand the choice xi = 0.08 was suggested comparing P{M) with 
the universal Ising curve, similarly to what has been done above for the Krypton- Xenon 
mixture). The simulation box linear dimensions were L = 9(jp, L = ll.Scp and L = IS.ScTp. 
For a fixed xi and fis, fJ'p is always fixed so that P{M) satisfies the equal area rule. Then, 
we compute second and fourth order cumulants B2 = {O'^) / -B4 = {O'^) / {O"^)^ (where 
O = AI — (M)) as a function of fig, for different L (fig. [T8k). It is seen that rather well- 
defined intersection points of the curves B2 and i?4 vs. fis for different choices of L do in fact 
occur at fis = —2.058 ± 0.001. Using a simple comparison with the universal critical curve 
(which is the method we have used for all the critical lines of fig. [T6l) we previously obtained 
a value fig = —2.06, which is very close to the value determined above (in particular within 
the 0.5% which is our general estimate of errorbars). The intersection for B2 occurs close 
to the theoretical valued = 1.22382, while the intersection point for the curves is 
somewhat too low. This may indicate that the choice of the mixing parameter xi is not 
optimal. However, varying B4 as a function of Xi, at fixed choices of fig we observe that i?4 
has a shallow minimum near the chosen value Xi = 0.08 (see inset of fig. [T8k). This could 
justify in some sense the apriori chosen value for xi and gives an estimate of the systematic 
error related to the choice of xi, which is very small. (We notice that variation of B4 with xi 
is not strong enough that one could tune the parameters such that the intersection occurs 
strictly at the theoretical value, and with a distribution function P{M) which is far from 
looking like the Ising curve). We conclude that all the sizes from L = 9ap up to L = 13.5(jp 
are not yet in the asymptotic region of finite size scaling, and so various corrections to finite 
size scaling occur which could only be disentangled if a much wider range of L were at our 
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disposal. In view of such possible systematic errors, we have allowed an error bar for the 
critical values of /i^ and yUp of 3 parts in a thousand, three times as large as one would 
conclude from a naive analysis of fig. [T8k . In any case, the uncertainties resulting from these 
finite size analysis are much less than the deviation between our model predictions and the 
experimental data. 
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Indeed using results of Ref. 



36l the CG parameters for given substances (with a reasonable 
quadrupolar moment) can be computed in a straightforward way without additional simula- 
tions. 

109 However sometimes is far from being clear if deviations from LB combining rules are infact 

ar substances 



42 were more 



compensation of some bias of the model used. For instance in the case of po 
the present work shows very clearly how strong violations of the rules in Ref. 
properly due to a bad modelling for the solvents. 
110 It is interesting to observe how atomistic models that are not improved to describe all the pure 
substances phase diagram but take as experimental input only the critical temperature and 
density (like in our case) are less accurate than our simpler (see for instance the discussion of 



EPM2 models in 



36|). 
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TABLE I: We report all the simulation parameters used in the present work and the critical 
parameters of the pure components that have been used. In brackets we also report the errors 
(if n digits are reported, the error applies to the same last digits). The errors for e and a have 
been estimated using the experimental errors for the critical points^^ and an error of 0.5% in the 
simulation critical points as estimated previously. As a consequence we note for instance that 
in the three models used for CeHg the three values for a are almost compatible with our error 
bar. It is important to observe that in this discussion we have disregarded the huge error in the 
quadrupolar moment Q (and as a consequence in q) . We refer to the text for a discussion of this 
point. 



FIG. 1: (Color online) Illustration of the coarse-graining procedure: in the case of hexadecane, 
three successive C-C bonds are integrated into one bead (dotted circle). The oligomer, containing 
50 atoms (or 16 "united atoms", CII3 or CII2, respectively) is thus reduced to an effective chain 
of 5 beads. Neighboring beads along a chain interact with a combination of Lennard Jones (LJ) 
and finitely extensible nonlinear elastic (FENE) potentials. Non-bonded beads only interact with 
a single LJ potential. Carbon dioxide is represented by a point particle, which carries a quadrupole 
moment. 

FIG. 2: (Color online) Coexistence curve for Kr (upper curve) and Xe (lower curve) in the temper- 
ature density-plane (a), and interface tension plotted vs. temperature (b). Broken curves indicate 
the experimental data,^^ asterisks MC results where pc^^ was used to adjust a, while circles show 
data where pf^"^ was used to adjust a (cf. text). 
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FIG. 3: (Color online) Coexistence pressures for pure Krypton (right) and Xenon (left) plotted 
vs. temperature. The broken line shows the experimental data for the pure noble gases, while the 
full curve is the projection of the critical line of the binary mixture (to be discussed in Sec. III.A).— 
The symbols denote our simulation data (the notation is the same as in fig. [2|). 

FIG. 4: (Color online) Coexistence curve for CH4 in the temperature-density plane (a), vapor 
pressure at coexistence (h) and surface tension plotted vs. temperature (c). Broken curves show 



experimental data (Ref. 
the MSA predictions. 



75I ). crosses the Monte Carlo simulation results, while the full lines show 



FIG. 5: (Color online) Coexistence curve describing vapor-liquid equilibrium for benzene (CgHg) 
in the temperature density plane (a), temperature dependence of the vapor pressure (b) and the 
interfacial tensions (c) at phase coexistence. The full curve is the result of fitting Tc^^, pf^^— to 
a simple LJ model without taking into account any contribution from quadrupolar interactions, 
while triangles are Monte Carlo results using Eq. Q and crosses Eq. (|lip . respectively. Broken 
curves are the experimental dataj^ 

FIG. 6: (Color online) Coexistence densities for the alkanes studied in the present paper. The fol- 
lowing substances are reported (from below): Methane (CII4), Propane (CaHg), Pentane (C5H12) 
and Hexadecane (C16H34). Curves are experimental data,— while the open circles are our simula- 
tion results. 

FIG. 7: (Color online) Coexistence pressures plotted vs. temperature, for the alkanes studied in the 
present paper, namely CH4, CsHg, C5H12 and C16H34 (from left to right). Curves are experimental 
data^^ dots show our simulation results. 

FIG. 8: (Color online) Interface tensions plotted vs. temperature, for the alkanes studied in the 
present paper, namely CH4, CsHg, C5H12 and C16H34 (from left to right). Curves are experimental 
data,^^-- dots show our simulation results. 

FIG. 9: (Color online) (a) Krypton concentration at criticality plotted vs. temperature; the curve 
shows the experimental data,— while symbols denote our simulation data (the notation is the same 
as in fig. E]). 
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FIG. 10: (Color online) Isothermal slice through the phase diagram of the mixture of CII4 and 
C3H8 at T = ?>TIK (a) and T = 271 K (b), using the molar fraction x of CH4 as abscissa variable 
and pressure p as the ordinate variable. Dots are experimental data^^ broken curves result from 
our TPTl-MSA approximation, and symbols denote Monte Carlo data (see text). Triangles are 
MC results for the critical points. 

FIG. 11: (Color online) Isothermal slice through the phase diagram of the mixture of CH4 and 
C5H12 at T = ?>1SK (a) and T = 237 K (b). Triangles denote experimental data,-^*^ broken 
curve denotes TPTl-MSA, and asterisks denote Monte Carlo results (which were only taken for 
T = 378K). Triangles are MC results for the critical points. 

FIG. 12: (Color online) Isothermal slices through the phase diagram of the CII4+CO2 system at 
T = 270K (a), T = 250K (b) and T = 230K (c). Dots represent experimental data^— while the 
broken curves are the results of TPTl-MSA when CO2 is represented as a point particle with 
no quadrupole moment {qc = 0). The full curves are results of TPTl-MSA with the spherically 
averaged quadrupolar interaction using the parameters of Eq. ([8]) {qc = 0.387). The triangle shows 
the MC result for the critical point. 

FIG. 13: (Color online) Isothermal slices through the phase diagram of the CII4 + CgHg systems 
at T = 501. 15A' (a), T = 461. 85K (b) and T = 421.05^ (c). Full dots show experimental data^ 
curves are calculations based on TPTl-MSA: broken curves denote the simple LJ model {qc = 0), 
dash-dotted curves are based on Eq. and full curves on Eq. (jlip . The triangle shows the MC 
result for the critical point. 

FIG. 14: (Color online) Isothermal slice through the phase diagram of the CO2 + C5II12 system 
at T = 423. 48i^ (a) and T = 344.34i(' (b). Full dots represent experimental data,^=2i asterisks our 
Monte Carlo results for the model, Eq. (jlOp . while the full curve is the corresponding TPTl-MSA 
prediction. The broken curve shows the corresponding TPTl-MSA result for a CO2 model with 
no quadrupole moment {qc = 0). Triangles are MC results for the critical points. 
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FIG. 15: (Color online) Isothermal slice through the phase diagram of the CO2+C16H34 system at 
T = 486-fC, showing MC results for the present model (open circles) and comparing them to the 
results of the previous simulations'^^ with = 0, ^ = 1 (full dots) and = 0, ^ = 0.886 (asterisks). 
Squares show two sets of experimental data^^ at two temperatures that bracket the temperature 
used in the simulation. Triangles are MC results for the critical point. 



FIG. 16: (Color online) Critical line of the CO2 + C16H34 mixture, projected onto the p*,T* plane 
(pressure p and temperature T are rescaled with the LJ parameters of the effective monomers of 
hexadecane as usual, p* = pe/a^ and T* = ksT/e). Different symbols (as indicated in the figure) 
denote data with gc = 0,^ = 0.886 (top curve) and ^ = 0.9,^ = 1 (lowest curve), as well as data 
for nonzero quadrupole moment, Qc = 0.387 and q = 0.47, respectively. 



FIG. 17: (Color online) Final normalized order parameter histogram P{M) of Xenon-Krypton 
mixtures (curve) at Ti = 228.78. The simulated (s) and reweighted (r) parameters (in units of 
and a^) are ^i]^^ = -2.254, = -2.2512, /i*^ = -3.972, ^u*^ = -3.9792, xi = 0.4, X2 = 0.03. 
The dots show the universal 3d Ising model distribution. 



FIG. 18: (Color online) (a) Plot of and B2 as a function of jig for the CO2 + C16II34 mixture 
(T=1.16). The chemical potential was always chosen such that the equal weight rule was 
obeyed. Three different box linear dimensions are included, as indicated. The broken horizontal 
lines indicate the universal values B\ and S| of the Ising model at criticality, where the intersections 
of the curves for ^4 and B2 in the finite size scaling limit (L ^ 00) should occur. These data 
have been generated for the mixing parameter xi = 0.09. The inset shows a plot of -B4 vs. xi for 
L = 13.5crp for three different values of /x^. (b) Probability distribution P{np,ns) of the numbers 
of polymers (up) and solvent molecules (n^) for L = 13.5(Tp at criticality. 
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